in 
o 
o 

o 

Q 

o 

(N 



Solution of the Poisson equation for two dimensional periodic structures (slabs) in an 

overlapping localized site density scheme 

F. Tasnadi^Q 

'IFW Dresden e.V., 
P.O. Box 270 116, D-01171 Dresden, Germany 
(Dated: February 2, 2008) 

Bertaut's equivalent electric density idea (E. F. Bertaut, Journal de Physique 39, 1331 (1978)) 
is applied to the case of two dimensional periodic continuous charge density distributions. The 
following derivation differs from what was introduced by Bertaut. The presented method solves 
the Poisson equation for the scheme of overlapping localized site densities with periodic boundary 
conditions in the (x,y) plane and with the general finite voltage boundary condition in the perpen- 
dicular z-direction. As usual the long-range potential is calculated in the Fourier space. For the 
Ku 7^ case a Fourier transformation helps to calculate the solution in a three dimensional periodic 
sense, while for Kii = the required charge neutrality is the starting point. For both cases suitable 
representations of the spherical harmonics are needed to arrive at expressions that are convenient 
for numerical implementation. In this localized density scheme an explicit relation can be derived 
between the finite voltage in z-direction and the ^-component of the dipole density. 

Keywords: slabs, Ewald summation, Poisson equation 
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I. INTRODUCTION 

The solution of the Poisson equation for periodic struc- 
tures has great practical importance in many areas of 
solid state physics, e.g. in band structure calculations 
or in molecular dynamics simulations. The Bertaut's 
' equivalent electric density' ideai, realized in the pseudo- 
charge or multipole compensation method developed by 
Weinerl£ is nowadays the standard approach in plane- 
wave implementations for three dimensional (3D) peri- 
odic problems. The Bertaut approach to sum up the 
long-range electrostatic interactions is based on the re- 
placement of the point-like multipoles at every lattice 
sites by non overlapping charge distributions with equiv- 
alent multipoles. In the Weinert method every site might 
have several multipoles and the crystal is divided up into 
atomic site domains f2 (muffin-tin sphere) and into an 
interstitial part. To calculate the electrostatic potential 
in the interstitial region (the long-range potential con- 
tribution) Weinert realizes Bertaut's idea by introducing 
pseudo-charge densities confined to the atomic volumes 
f2. Inside £1 the Dirichlet problem for the original charge 
density is solved, using the potential value at the bound- 
ary of f2, which result from the interstitial solution. 

Dealing with 3D crystalline extended systems the pe- 
riodic Born-von Karman boundary conditions (PBC) 
in each direction are mandatory. For films or single 
slab geometries^ the PBC is applicable only to the ex- 
tended (x, y) (in-plane) directions. In the perpendicular 
z-direction the system has only finite extension and re- 
quires a different boundary condition. The most general 
physical situation allows for a finite voltage at infinity 
in z-direction, which would be the case for slabs without 
a in-plane mirror symmetry. This general finite voltage 
condition in the z-direction coincides with the physical 
picture of a charged capacitor. 

A possible treatment of films is the super-cell (vacuum- 



film-vacuum) method with 3D periodicity, which de- 
scribes the electrostatic potential inside the super-cell as 
a continuous function and thus excludes any step like 
discontinuities. The introducing of an artificial dipole 
layer in the vacuum region^ allows to model the finite 
voltage situation. However, there is still a possibility 
of spurious image interactions. Furthermore, an addi- 
tional parameter, the vacuum thickness, has to be con- 
verged. The plane-wave based layer-FLAPW method 6,7 
treats the single slab geometry as a real 2D periodic unit 
cell extending from — oo to +oo in z-direction with three 
different regions: the muffin-tin, the interstitial, and a 
vacuum (no charge) segment. The method is applicable 
to surface problems without any invocation of an artificial 
periodicity. In this method a finite voltage may occur be- 
tween the two interstitial boundaries. It seems that this 
treatment of the potential step feature via plane-waves 
produces thicker interstitial regions than necessary. 

This work presents the general solution of the Pois- 
son equation for two dimensional periodic systems based 
on the multipole compensation method for a scheme of 
overlapping localized site densities. In the compact site 
density ansatz (see Refs&S) the total electronic charge 
density is written as a locally finite sum without any in- 
terstitial density. The presented derivation differs from 
what is given by Weinert or Bertauli for 3D periodic 
systems. The correct capacitor-like boundary condition 
is taken explicitly into account, which results in a direct 
relation between the finite voltage and the dipole of the 
slab as is expected from macroscopic electrostatics of a 
film geometry. 

In Appendix iDl we shortly sketch how our ideas apply 
to the ID periodic problem, which differs in methodology 
from the treatmentsiSiii*^ presented elsewhere, however 
resulting in the same expressions. 
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II. THE MULTIPOLE COMPENSATION IN 
THE LOCALIZED SITE DENSITY SCHEME 

For crystalline structures in the thermodynamic limit 
there is an infinite periodic adjustment of the finite num- 
ber of sites in the unit cell {si, . . . , s<j}, by the Bravais 
vectors R. In the case of film geometry the Bravais vec- 
tors R span one of the 2D Bravais lattices. In other 
words, the slab geometry means a three dimensional ob- 
ject with two dimensional periodicity. The symmetry 
groups of this geometry are called layer groups 3 . Ac- 
cording to the main assumption at any point r in the 
extended solid only a finite number of site charges with 
compact support f2 s . j £ {1, . . . , d} will overlap, as it is 
shown in Fig. ^ Each site domain may contain several 
neighboring unit cells. 




FIG. 1: Schematic site charge (compact support) distribution 
are overlapping at the point r. 

Using the local expansion ansatz 

m^^lfo/lc^e*^' (1) 
R 3 

for the Bloch state |kn), where k is the crystal momen- 
tum, n is the band index and C is a composite index for 
atomic orbitals containing the principal and the angular 
quantum numbers, a local expression can be derived for 
the total density p 

occ d 

p(r) = 5>|kn)(kn|r) = ££>(r - R - s,), (2) 

k,n R 3=1 

which defines the site-density pj and the corresponding 
domain f2 Sj +R, see Koepernik and Eschrig^. The total 
charge density 

rGO R+Bj 

P(r)= £ ft( r_R_s i) = 

rGO R+Sj 

]T J2p jL (\r-R- Sj \)Y L (r-K- Sj ), (3) 

Rj L 



is accordingly a finite sum of site contributions and is 
written as a uniformly convergent sum of spherical har- 
monics Yl, L = (£,m). The first line defines the total 
density by the site densities. The complicated restriction 
in the summation allows only such R and j indices that 
the domain of the corresponding site density contains 
the point r. The lattice periodicity of the total electronic 
density p(r) is not broken because of the r dependence 
of the summation. 

To calculate the electrostatic potential 4>(r) of the sys- 
tem one needs to solve the Poisson equation 



A 0(r) 



ref2R, +Sj . 

£ 

R,i 



p j (r--R-s j )-y i Z j 5(r-R- Sj ) (4) 



R,j 



at any r £ R 3 , where no restriction is used in the second 
Dirac delta summation. To get 'bulk' quantities PBC is 
used in the periodic directions (say x and y) while in 
the third direction the general finite voltage boundary 
condition is considered. Trough the paper atomic units 
are used. 

We now introduce the generalized Ewald densities ex- 
panded by spherical harmonics 



p EwaM (r)= £ p E„ald (r)yi(r) 
L 



PiL {r)=Aj L Mere ' , Ale 



2p 



,2£+3 



Ewald 



L 



T(l+I) 
(5) 



where p is the usual Ewald parameter and the multipole 
freedom is settled in Ajl . Afe is the normalization to have 
unity multipole components for the Gaussian without 
AjL. With that normalization Aji, are not dimension- 
less quantities, but their dimensions are [Ajl] = metcr f . 
The last line of Eq. JSJ) indicates the most reasonable cou- 
pling of terms what is used in the calculations below. By 
the additional multipole compensators (Ewald charges) 
a modified electrostatic problem is created with a new 
Poisson equation 

A 4>(r) = p(r) 

p(r)= ]T p i (r-R- Sj )-£^(r-R-s 3 -)+ 



Rj 



Rj 



^pf^V-R-s,), (6) 
Rj 

at any r £ R 3 . Taking the difference of the two Pois- 
son equations and using the linearity of the differential 
operator A, the Ewald problem 

A Ewald (r) = pf WaM ( r - R - s j ) = P EwaU ( r ) . ( 7 ) 



3 



is derived. Thus, the complicated original Poisson prob- 
lem is subdivided into two simpler problems, 4>(r) — 
4>(r) — ^> Ewald (r). Before discussing the boundary con- 
ditions for the latter two equations the solutions of the 
multipole compensations are needed because they define 
the Ewald densities and thus the two Poisson equations 
Eas. l|6l7|) . Keeping the 2D periodicity in both subprob- 
lems is compatible with the original PBC. 

A difficulty shows up in the solution of Eq.© by the 
non-restricted sum of overlapping non-compact Ewald 
densities. To get rid of this difficulty one can assume 
that the Ewald densities are negligible outside the site 
domains. This results the same finite number of terms in 
the last sum of Eq.|@ as in the first sum. Using the ap- 
proximate compactness of the Ewald densities the multi- 
pole compensation equations are written as integrals over 
the whole space 







/ d 3 rY; m (ry(p J (r-R-s J )-Z ] 6(r-K-s,) + 

(r-R-sA (8) 



Ewald / 

Pj 



The introduced assumption simplifies the Poisson equa- 
tion Eq.© into the form 



A </>(r) w p(r), 



(9) 



where the r.h.s. is given as an r-dependent finite sum 
of site densities. This equation is referred below as the 
short-range Poisson equation. The accuracy of the ap- 
proximation is measured by the localization of the ap- 
plied Ewald densities. 

The linearity of A allows one to look for the solution 
also in a finite sum 

4>(v) = ^'(r-R-Sj), (10) 

Rj 

at any point r, where 4>j has the same compact support 
as the site density pj. Accordingly on each site domain 
f2 Sj a Dirichlct boundary problem is derived by the multi- 
pole compensations. Extrapolating these boundary con- 
ditions to the discussed equation Eq.@) gives that <j) has 
to vanish at the boundary of a composite compact region. 
Of course, the in-plane periodicity is kept by means of the 
r-dependent summation. The solution of a site boundary 
problem with vanishing potential at the boundaries is not 
discussed in the paper because it corresponds with the 3D 
periodic case which is already published in Refs^. 

Accordingly, the Ewald problem shows also the obvi- 
ous infinite lattice periodicity and carries the finite volt- 
age boundary condition in the surface normal direction. 
The infinite lattice periodicity allows one to use Fourier 
expansion for the functions 

p Ewald (r) = ^e lK n r p Ewald (K||,z) 

K|| 

i r Ewald (K| h z) (11) 



iEwald 



(r)=J>«- 



and the Ewald problem Eq.lJTJ) is transformed into the 
reciprocal K|| space 



A_0 E » ald (K||,z) -40 Ewald (K| h z) = p Ewald (K| h z) 
p Ewald (K||,z) = l^ / d%T pf™ W (r-Sj>- iKuT , 



(12) 



where U means the volume of the 2D unit cell and d 2 r 
denotes dxdy. The structure of the last equation dictates 
different treatment for Ku = and K\\ ^ 0. Since the 
A'n =0 Fourier coefficient is a constant in- plane function 
at any z value the A| | = equation provides the required 
finite voltage boundary condition whereas all the K\ | ^ 
solutions vanish at z = ±oo. 



III. SOLUTION OF THE EWALD PROBLEM 

The Ewald boundary problem is defined in the previous 
section with two subproblems. As mentioned previously 
the A|| = case differs notably from the ordinary K\\ ^ 
case and thus needs a different treatment. The K\\ ^ 
case contains a non-singular differential operator on the 
l.h.s. which ensures applicability of the Green function 
technique in the solution. On the other hand the K\\ = 
case can be solved by direct integration starting from the 
charge neutrality condition and using the given boundary 
conditions. 



A. Solution of Eq.JT^ for K N / 

This equation is a so called ID modified Helmholtz 
equation and its Green function reads 



G(z,z) = — 



e -\z-z\K n 



2 An 



(13) 



including the vanishing boundary conditions at ±oo, 
which results a complicated integral for the solution 



/Ewald 



(K,,,«) 



e -iK|| Sj - 



3 



u 



,-\z-s jz -z\K\ 



2Ku 



^Ewald /- 



(f) e 



-iKnf 



(14) 



K 



The difficulty arises with the spherical harmonics expan- 
sion of p Ewald , which is the most reasonable choice to 
have a straightforward solution of the short-range Pois- 
son problem. Here this advantage turns into a trouble- 
some disadvantage because the Green function and the 
exponential e~ lK n r are preferably handled in a Cartesian 
coordinate system. 
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With the help of the Fourier transform of the Green 
function corresponding to the z — Sj Z variable , 



1 



n/2tt J- 



+ 00 



F(cj,z) = 



-\z-sjz-z\K u 



2Ku 



2n{u: 2 + K 2 ) 



(15) 



one can introduce an e~ tuJZ function which can be coupled 
to the e _zK n r function. The introduction of the new three 
dimensional vector 



KU-, = («II 



(16) 



allows one to make use of the Rayleigh relation 
e -«(o0? = 47r^(-0'Vv(K(a;))j>(^(^)^)^(f) 5 (17) 



where ji are the spherical Bessel functions. In this way 
the problem is transformed into a form that is similar to 
the 3D periodic case with the only difference that now 
the third reciprocal direction u> is continuous. 

Interchanging the order of the Fourier J tko and the 
real space integration J d 3 r and doing a straightforward 
but lengthy calculation the following equation is resulted 



/Ewald 



(K||^) = -£- 
K, 



,-iK, 



(-*)< 



u 



e ^AjLp x 



duj 



(2 P y 



(18) 



Here, the prefactor in the first brackets is dimensionlcss 
and the w-integral shows the dimension of length. In the 
calculation of the w-integral a representation of the spher- 
ical harmonics by hypergeometric functions is used. An 
appropriate choice provides not only a simplified writing 
of the solution but also a representation that allows an 
accurate numerical implementation. It is worth to men- 
tion that taking any representation of Yj, in Eq. (|18|) the 
L = (0,0) case will result the same integral, this case 
is come down from the Fourier transform of the Green 
function. According to the handy relations between the 
spherical angle-coordinates cp) = f and the variables 
K|| and u> 



sini9 



(K\\ x + iK ny y 
KIT 1 



cottf 



(19) 



the most reasonable choice is given by 
2 e+1 ^Y L (r) C(t,m) 



(2£+l)\l 
{ 2^1 ( 



Yjn t-rn . 21-1 
2 ' 2 ' 2 ' sin 



e^sintf)* x 
1 <i 



til- 



m is even 



cot 1? F(- 



l+m-l e-m-1 . 21-1 . 



1 <| 
7~d> 



(20) 



2 ' 2 ' 2 ' sin- 1 -0 

if £ + m is odd 



(see Varshalovichi^ , page 137, Eq.(30)) with the real 
valued quantity 



C(£,m) 



2*(-l)L^J+i 



y/(l + m)\(l-m)\(2e+l)' 



(21) 



where [• ■ • J indicates the greatest integer function or 
floor. 2-F1 is the hypergeometric function and (•••)■■ de- 
notes the double factorial defined by 




2) ... 5 • 3 • 1 n > , odd 
2) ... 6 • 4 • 2 n < , even 
1 n = -1,0 



(22) 



In Ea. l|18|l the Ajl quantities arc calculated by means 
of the multipole compensation. Experience shows that 
the quantity Ajl decays quickly with higher £ values. 
This fact enables one to apply a cut-off in the summation 
which usually has a value of £ max = 12. 

For any finite I the hypergeometric functions above 
terminate after finite number of terms 



c/o 



iax (£,m) 

E 

k=0 



e/o / 



a k (£,m)K- 2k (K* 



,2\k 



(23) 



because one of the first two arguments of F is always 
a non-positive integer, k > 0. Here the subscript e/o 
distinguishes the even and odd values of (£ + to). The 
expansion coefficients o^/° (£, to) are real with alternating 
sign and /c ma x(^, m ) is the £ and to dependent termination 
value, see Table [|] 



(l,m) 




(i,m) 




(l,m) 


femax 


(£,m) 




(0,0) 





(3, -3) 





(4, -2) 


1 


(5,-3) 


l 


(1,-1) 





(3,-2) 





(4,-1) 


1 


(5,-2) 


l 


(1,0) 





(3,-1) 


1 


(4,0) 


2 


(5,-1) 


2 


(1,1) 





(3,0) 


1 


(4,1) 


1 


(5,0) 


2 


(2,-2) 





(3,1) 


1 


(4,2) 


1 


(5,1) 


2 


(2,-1) 





(3,2) 





(4,3) 





(5,2) 


1 


(2,0) 


1 


(3,3) 





(4,4) 





(5,3) 


1 


(2,1) 





(4, -4) 





(5,-5) 





(5,4) 





(2,2) 





(4, -3) 





(5, -4) 





(5,5) 






TABLE I: The k m ^(£,m) values for £ g {1, . . . , 5}. 



Henceforth the even and odd contributions are treated 
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separately. The w-integrals from Eq. l|18|l with the third 
bracketed prefactor are written as 

T =C{i,m)e^P tm K\\ \ dio 1 e »>b-i*), 

(K,, + LO ) 



odd. Taking the values of fc max from Tabled it turns out 
that in the cases of fc max (£, m) = only X 3 and Zj provide 
contributions. These two integrals are defined as 



R 



1° = C(£,m)e mip P em / duj uj 



23(«,i,-K"||) = / dx - 

JR 1 



-if 2 z 2 /(4p 2 ) 



,iX||(z-Sj z )a 



Zi(z,j,Ku) = / da; 



2p V 2 P 



xe 



+ x' 



iKu(z-Sjz)x 



(24) 



R 



1 + X 2 



(28) 



Inserting Eq.J^SJ into the above integrals results the fi- 
nite series 

2i? e/° e -a;7(4p 2 ) g-^ 2 /(4p 2 ) 



(if 2 +^) (iffj+a, 2 ) 

^|7 2 ~ 2 '^ 2n e-" 2/(4p2) ^ o/o (n,£,m) (25) 



in the integrands, where the n-summation comes from 

Ml 



the binomial expansion of (K?, +Lu 2 ) k with n > and 



fcmax(<,m)-l 



g c />,^,m) = Yl J" ( 26 ) 

To derive Ea. 1)25(1 the order of the finite fc and n sum- 
mations is interchanged. After substituting the last two 
equations into Ea. l|24|l a convenient handling, 

T = C(£,m)e lmv P tm x 

fema*(<,m) — 1 

X 3 (z,j,X||) + £ Z 2n e? e (n,i) 

n=0 

1° = C(£,m)e lmip P em x 
2k(z,i,lf||) + £ J 2n+1 g°(n,L). (27) 

is achieved for the lengthy expressions by introducing 
three integrals 13,24 and 2" 7 , where 7 can be even or 



where u = K\ 1 x is used. They look like a Fourier transfor- 
mation which suggests the use of the convolution theorem 
of Fourier transformations. The final results of the inte- 
grals are presented in Appendix [S] The T<z n and Tm+i 
integrals are treated together by introducing the dimen- 
sionless new variable y 

2 7+1 / dy y~ i e~ {K » y2/p2)+2{lK u [z ~ Slz))y = 



,•7 



2p 



7+1 LtJ 



L IU t=0 



(29) 



where u> = 2K\\y and £ — 2 > 7. The factor in the first 
bracket describes the if 1 1 dependence and the exponential 
from Pi m ensures its decaying behavior for large Kit. It 
is also worth to mention that the parity of 7 completely 
separates the results. For even 7 the integral has a real 
value while for the odd cases it is imaginary. 

From the numerical numerical point of view in the 
above sums the the quickly vanishing Aj^s determine the 
order of the contributions and ensure the inclusion of all 
relevant contributions. 

A partial simplification can be achieved in the T c l° 
calculation if the order of the n- and t-summation is in- 
terchanged. Hereby the n-summation with the [z — Sj z ) 
functions goes ahead and the same upper limits can be 
used in the summations, 



I c = C(l,m)e im *Pz m l 3 + 

C(l,m)e im *P em J2 {-l) n e- p2{ *- 8i ° )2 W-s jz )? n f (jf) (Y?*W - n + h) G c (t, L); 

n=0 t=n \ II ' V V-"/ 

1° = C(£, m)e mv P lm l A + 



iC(£, my^P tm ^(-iJ-e-^C-i-)' (p{z - s JZ )f^ ""ff (jf) ^ ( r(* - n + 5) ) G°(t, L) 



(30) 
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The final real valued ifii ^ Ewald potential can be 
easily derived from the result 



lE r ld (r) = EE'E( e ' K| 



o if II 



C(£, to) 



Ira 



Z 3 + En J 2nG C 






to even 



£ + to odd 
(31) 



where the two dimensional column vector denotes the 
real and imaginary components of an arbitrary complex 
number. 

Simpler expressions can be found for the special cases 
if (z — Sj z ) = 0. It is easily established by symmetry 
reasons that only for even 7 gives the integral I 7 (see 
Ea. (|29H ) a non-zero result 



(32) 



For the other two integrals (23,1.4) the results can be 
calculated by using the general expressions given in the 
Appendix ^] (see also Gradshteyni4, 3.466/1.) 



M z = s jz ) = 7re*"n /(4p2) Erfc 



^1 
2p 



l±(z 



0, 



(33) 



which tells that only the even (i + to) terms contribute. 

Easy to check that all the integral results vanish at the 
boundaries z = ±00 which ensures that the final formula 
fulfills the required boundary conditions. 



B. Solution of Eq.JTJJl for K = 

In this section the solution of the boundary problem 

-^ Ewald (0,z) = p Ewald (0,z), 
dz z 

^Ewald^ +QO ) _ ^Ewald^ _ QQ ^ = ^ (34) 

is looked for, where the Ewald densities are defined 
through the multipole compensation and C2 refers to 
the already introduced finite voltage. Taking only the 
monopole equation 

/ dh ( Pj (r - Sj - R) + p Ewald (r - Sj R)) -Z s = 

(35) 

together with the charge neutrality equation 

/+00 p 
dz I dxdy Pj {v - Sj - R) - ^ Zj = (36) 
it,, -°° * u 3 



the well-known fact of the neutrality of the Ewald density 
lattice can be derived 

£ / d?v pf^(r - Bj ) = «/X £ A m°) = °- (37) 
3 Jn V 3 

For a 2D periodic system Ea. (|37ll has the consequence 
dz ap Ewald (if 1 1 = 0, z) = 0, a E R \ {0} (38) 

which can accurately be approximated by a bounded 
integration range [D\ OWl D up ] for well localized Ewald 
charges. In words, one uses the bounded region where 
the Ewald charge distribution is non-zero. Integrating 
Eq. (|34|) and applying the previous equation results 



dz 2 



dza[^ Ewald (O,2 



0. 



(39) 



which yields the boundary behavior of the negative elec- 
tric field component E(z) = d0 Ewald (O, z)/dz 



dz 



dE(z) 
dz 



a(E(D up )-E(D low ))*0. (40) 



Applying the Gauss law outside of this bounded region 
results the boundary problem for E 

dE(z) 



dz 



= P 



^Ewald 



(0,z), 



E(z < Aow) = const. 
E(z > D up ) = const., 



(41) 



which defines another first order boundary problem for 
the if 1 1 = Ewald potential 



d0 bwald (O,z) 
dz 



= E{z), 



^Ewald 
/.Ewald 



(0, z < Aow) = 
(0, z > D up ) = C 2 



, (42) 



The finite values at the potential boundaries dictate the 
const. =0 conditions for E. It is worth to mention that 
applying non-zero constant boundary conditions for E 
contradicts with the constant behavior of the potential 
outside of the bounded range [Aow, D up }. 

Both problems can be solved by simple integration 
which results the final expression 

Ewald (O,z) = 



dz' E{z') 



D u , 



dz' E(z')\+^ 



Aow < z < D up (43) 
with explicit inclusion of the potential difference C2. 

C. Evaluation of Eq.ll43ll 

Integrating the boundary problem in Eq. l|41[) results 
E(z') = 

f dx"dy"p*r ld {r")Y L {v") (44) 
Jut? 



1 E 1 



dz" - / dz" X 

D Jz'-su I 
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with a common D, for which D — > oo, and E will satisfy 
the requirement given by Eq. lJ58|) . Unlike in the K\\ ^ 
case the above integral allows one to use cylindrical coor- 
dinates (g, ip, z"). Using the representation of the spher- 
ical harmonics (see Varshalovich 1 -, page 137, Eq.(33)), 



„, I- \m\ I- \m\ - 1 , , , „. 

XF( J-l, LJ ; | m | + 1; _ tan 2 0) 



2 

CmO 



|to|!2I™I 

- tan 2 ■& 

(-l) m m > 
1 m < 



(45) 



immediately results 2nS m fi for the (/^-integral and simpli- 
fies the expression for E to 



Here, the explicit form of the hypergeometric coefficients 
a is applied and C ( j, £, k) is defined as a dimensionless 
quantity. The one dimensional integral has the dimension 
of a surface charge [qt,j] = 1/m 2 . In the expression above 
the [• • • J bracketing is used again to denote the floor 
function. 

If t is odd (t = 2A + 1) then 



<?2A+lj 



n=0 



(48) 



(see GradshteynM, 3.351/2.). For z' — Sj Z only the n = 
term is non-zero. For even t (t = 2A) the result 



it 



LiJ 



dz" z "«-2*0 e -*"V X 



2fc+1 - e 2 



dg g ZK+l e 



(46) 



A = g ,j 



?2A,j = < 



Erf((z' - s jz )p) 
A > gsxj = f(\)q 0tj - ^e-(*'- s -> 2 *> 

Eti (/(A)//W) ((*'-^)p) 



x 

2ra-l 



(49) 

is proven in Appendix iBl and expressed with the help of 
the error function Erf ( • • • ) and a gamma function 



Here the hypergeometric function is written again by a 
finite sum with coefficients &k{l). After some algebraic 
manipulations with finite I summation E is expressed as 



A = 



(2A-1)!!/2 A A>0 



(50) 



^( Z ')=E E E C(j,2X + 2k,k) 

j=l A=0 fc=0 



^m,x-(2A+l) j 



E E ^A+u E C(i I 2A + l + 2fc,fe) 

j=l A=0 fc=0 

~+D 



1 



dz - 

— D J z' —Sj 



dz P 3 (z P fe- z p 



C(j,£,k) 



p 2 w V2TTTV 2 Vr( 



■ff. (47) 



which shows also the identity 



£/(' 

n=0 



A! 



2/(A+l) 



(51) 



proven by induction to A. This identity is used below in 
the calculations of the boundary properties of E. Eq. l|4^|) 
is also valid in the case of z' = Sj z . Further algebraic 
manipulations on E leads to the expression 



E(z') = ^ Erf ((^' " s jz )p)P (j,p, 0) - £ 

3=1 3 = 1 



2 „ 2 a L — 2 — J 



X] Pe(j,P,n)((z' - s 3Z ) P f n 



n=0 



E 

3=1 



P - ]T P (i,p,n)((z'- Si ^) 2 "- 1 , 



(52) 
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where 



2 L ^ 1J 2A M2k+m (p/2) 2k + 1 y7 r(2fc + 2) / * (-l) k - x (2A)! 



)_ nl ^ P 2 W v^4fcT3rf2fc+ll U- 



P 2 W VifcT3r(2fc+|) ^ fc! (2A + 1)!!VA 

p r ,- „ = J_ L V J 2^-, (2M) ( P /2) 2fe ^ r(2fc + i) / * (-i) fc - A (k ' 

j /(n) ^ P 2 W V4fc+Tr(2fc + I) ^ fc! 

2A 

Po(j,P,0)=C(j > J 0) = -^>. (53) 



The A; sums above are numerically feasible representations. Summing up P o (j,p,0) from Ea. i|53|) to all sites j 
results zero by the Ewald neutrality Ea. (|37|) . This identity is used to check the boundary properties of E. A detailed 
justification of the vanishing site sum of P (j,p, 0) starting from Ea. (|47|l can be found in Appendix IO The whole 
space integral of E connects the boundary values of Ewald (O, z), 



/+oo d . — 2 r+D 

dz'E{z')= hm V^^P o (j,P,0) / Erf(( 2 ' - s jz )p)dz' 
-oo D ^+°°~l 2 J-D 



2=1 

E Ye^/wV E J C'(i ) 2i + l li -A) = -V^^ S ,^, (0>0) -^^C'(i,l,0) 



j=l A=0 \n=0 



2=1 



2=1 



47T x I 4tt \ -k 

2^ S J*4?(0,0) ~ 77 V T A 2(i.o)- 



2=1 



2=1 



(54) 



and results the value of C2. For the third sum in the 
second line one can apply the identity from Ea. H51l) and 
the equation Ea. (|C6fl from Appendix IO 
The equation 



expresses the finite voltage C2 as the z component slab 
dipole (surface, area) density and gives the expected 
macroscopic electrostatics of the system. 



/ d 3 r zp.j(r — Sj) — dz dxdy zp{v) (55) 
j=i "'k 3 Jn Ju 



derived by the fact that the Bravais vectors have zero 
z components can be used in the z component dipole 
compensation. The resulted relation 



dz / dxdy zp{v) 
Ju 



/4tt x - 



E s i zZ i 



2 = 1 



2 = 1 



^E s J' z A>(o.o) = C 1 U + E 8 * zZ i 



(56) 
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Turning to the calculation of the potential by Ea. i|43[) one arrives to the final expression 



*Ewald 



(0, Z) = %PaU,P, 0)V5F((« ~ «i.)p)Erf((2 ~ S 3*)P) + E f P °0',P> 0)( 



3=1 



3=1 



E E (/WeC7,P,n))|V^lrf((«- ai ,)p) 



« L 2 J n // \ \2f— 1 

£ £ ('Wf.KM))K M ' ,v E fc f L - 



3=1 « =1 



E E ((n-l)lPo(j,P,n))- 



3=1 "=1 



^ t\ 2 

t=0 



(57) 



A straightforward calculation leads to the observation 
that only the first term and the third one in the second 
line have contribution at the boundaries. The calculation 
of the boundary values yields the same results what is 
already calculated in Ea. (l54ll . 



APPENDIX A: CALCULATION OF THE J 3 , T 4 
INTEGRALS 



By the Fourier transformations 
/ dx — 1 



2tt 7,r 1 + 



irx I " — It I 

e = ■» / — e 1 1 



SUMMARY 



The paper proposes a general approach to solve the 
Poisson equation for two dimensional periodic systems. 
The method is based on the pseudo-charge method^ and 
uses a localized site density scheme&S, where the total 
charge density is written as a locally finite sum. The 
bounded domains f2 Sj coincide with the compact support 
of the local electronic site charge densities pj(r — Sj). Go- 
ing beyond the Ewald technique, the approach not only 
compensates the site monopoles but also all higher mul- 
tipoles by the additional generalized Ewald site densi- 
ties. Solutions of the multipole compensation equations 
yield the multipole freedoms and unambiguously define 
the Ewald densities of every site. 

The proper boundary condition for the slab geometry, 
i.e. the finite voltage at the infinite boundaries in the z- 
direction, results in the correct macroscopic electrostat- 
ics. The finite potential difference at the boundaries is 
expressed by the total dipole (surface, area) density. To 
our knowledge, the explicit treatment of the general slab 
boundary conditions presents the most simple and direct 
way for self-consistent calculations of the finite voltage. 
The expressions presented here, may easily be applied to 
the electrostatics of lattices with point-like multipoles for 
molecular dynamical simulations. 

In Appendix [D] we shortly sketch how our ideas apply 
to the ID periodic problem, which differs in methodology 
from the treatmentsiSiii*^ presented elsewhere, however 
resulting in the same expressions. 



/2vr 



dx e V ^ J ' e ^- T > 



r)x _ PV2 — 



TR 



1 

e il 



1 

/2-7T JJR, 



dx xe V 2p J " e l( - s - 



_ T)x = iAp 3 (S-T) i 



(6—r)*p-° 



(Al) 

(see Gradshteynii 3.354/5., 3.323/2., 3.462/2.) ,where 
5j = K\\(z — Sj z ) and by the convolution theorem one 
can calculate the required integrals 

' f drT[f](r)T[cM5 3 ~r)=T[fg](5 j ) 



V27T JJH 

71" ft- 2 /tA„ 2 



|e^i/(V) ( ^ (5 . )+ ^ ( _ 5 . )) 



where 



V{8) = e^Erfc 



K \\ P$ 
2p K\\ 



(A2) 



(A3) 



APPENDIX B: CALCULATION OF THE q 2X 
INTEGRAL 



Easy to find out the recursive relation 



<?2A, 



2p 2 



2A - 1 



2p 2 



"92A-2, 



(Bl) 



10 



for the needed integral. Using this relation one can fore- 
bode the final result presented in Ea. l49|) and prove it 
easily by induction to A. 



APPENDIX C: SPECIAL SUMS FOR THE 
BOUNDARY BEHAVIOR 

The P Q (j,p, 0) quantity is the coefficient of the error 
function resulted from the integral q2\,j, see Eq.(??). By 
this observation one can do the calculation on the follow- 
ing way: 



Checking the boundary behavior of the Ewald poten- 
tial in Ea. (|54|l one needs the equation 



L 2 \ t 



E E ^ A + 2f + l,i — A) = C(j, 1, 0)/(l), 



t=o x=o 



(C6) 

which can be proven on a similar way given above using 
the 



n!=ifo+ §-*)(*+i-o 

(t+f)(i + l)(i-A)!(i-A)! 



a t -x{2t+l) = 



(C7) 



EE E C(j,2X + 2k,k)f(X) 

j=l X=0 k=Q 



d [~ 2~ J t 



EC(i, 0,0)+E E J2C(j,2t,t-A)^± 



j=l t=l x=o 



2A 



J (00) 



p 2 U 



= E 

3=1 

j=l t=l r ^ A=0 

(CI) 



where 



S(t,X) = f(\)a t -x(2t)(t - A)!(-l) ; 



a t -A(2t) = 



mzo(t+i-i)(t+^-i) 

(t + i)(t + h(t-x)\(t-xy: 



(C2) 



The first term in Eq. (|Clfl vanishes by the charge neutral- 
ity conditions stated in Eq. (|37|l . Using the above expres- 
sion for a t -x(2t) one gets 



A = 



s(M) = ir|?^§(-i) A a>o 

and then one can easily establish that 



(C3) 



t is odd: S(t, A) = —S(t,t— A) A = 0,...,*=± 
i is even: 5(4, A) = i - A) A = 0, . . . , ~ — 1. 

By this property for odd t it is trivial that the last sum- 
mation in Ea. (|Cl|l gives zero. For even t (t = 28) 

X-S(26X) {4S - 1)U |f W (4J ' 1)!! ( 1)* 



A=0 



(46 - 1 )!! 
4^ 



A=l 



E( 2 A 5 )(-l) A -0. (C5) 

A=0 ^ ' 



expression for the expansion coefficients of the hyperge- 
ometric function in Ea. (|45|l . 



APPENDIX D: THE ONE DIMENSIONAL 
PERIODIC CASE 

This appendix shortly sketches how the method inher- 
ited from the above presented treatment apply to the one 
dimensional periodic case. The general solution, char- 
acterized by the multipole moments of the unitccll and 
based on the Poisson summation formula, is already pre- 
sented in the literature see Refsii2*ii. The Ewald method 
for point-like Coulomb and dipole-dipole interactions was 
worked out first by Portoi^. 

Fourier expansions in the Poisson equation yields the 
2D modified Helmholtz equation 



8P_ 

dx 2 



8P_ 

dy 2 



lEwald 



. . , (D1) 

The equation requires different treatments for K ^ and 
for the case of K — 0. The K ^ case is solved with the 
help of the Green function 



G(K,p-p') 



1 

2^ 



K (K\p-p'\), 



(D2) 



where Kg is the modified Bessel function of the second 
kind (K n ) for n — 0. Unlike in the 2D periodic case now 
one can use the Laplace transform of the Green function 



dx' dy' 



1 



TR 2 
d 

E 



dt- 



o+ 



L2tt / 

1 

2t 



-t(x-x') 2 -t(y-y') 2 



R 



dz' pf wM (r' - sj)e 



%Kz' 



(D3) 



and then apply the Fourier transforms (t > by the 
Laplace transform) 



-f(r-r') 2 



/2tt 



-u: 2 /(4t) 



t iui(r— r') 



R 



(D4) 
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for the two exponential terms in the bracket. Of course, 
the dt integration behind the uj x and u> y integrals can 
also be executed and results the two dimensional version 
of Ea. l|18[l . but its calculation turns to difficulties after 
applying the already introduced multipole shaped Ewald 
densities. Executing the integrations in a different order 



To obtain the solution of the equation with K = one 
can follow the method of separation of variables after 
some art of manipulations. The equation is rewritten in 
cylindrical coordinate system and the solution, according 
to the assumption of the finite number of non-zero elec- 
trostatic multipole moments, is looked for in the form of 



dt 



duj x 



R 



duj„ 



R 



'■'>.,' 



R ! 



(D5) 



allows one to introduce the Fourier vector K T {uj x , aj y ) — 
(uj x , ui y , K) similar to that is used above in Ea. H16(l . The 
cylindrical arrangement of the system suggests the appli- 
cation of Eq. (|45() to represent the spherical harmonics in 
the calculations. The phase factor in Y^ m can be given 
as 



(ui x + isgn(m)m y y m \ 

n\ m \ 



(D6) 



which leads after its binomial expansion to the general 
result of the duj x diOy integral 



l af3 {t) = ( z )«±2mod( Q ,2) 
^•\/3±2mod(/3,2) 

where 



d" 



d{x-s jx ) a 



T x {t) 



d(y- 



ly(t) , (D7) 



: exp 



p 2 t(i - sji) 



t 



with its special case 



i e {x,y} 

(D8) 

(D9) 



Having these integrals one should turn to the dt integra- 
tion, which shows the following form 



oo exp 
dt 



/ K 2 K 2 
\ it Ap 2 



At 2 



(D1Q) 



In general the dt integration leads to the so called leaky 
aquifer function W n . Especially for the case of L = (0, 0) 
- only monopole compensation - where a = [3 = the 
zeroth order function Wq is obtained 



dt 



4p 2 



np 2 t 



At 2 



t + p 2 



exp 



(-^[^(--^) 2 +P 2 (y-^) 2 ]) = 



-Wo {K 2 /(4p 2 ),p 2 (x ~ Sjx ) 2 + P 2 (y - s jy ) 2 )(Dll) 

Appropriate numerical representation of the leaky aquifer 
functions is discussed by the F.E. Harrisi£*i& and E.S. 
Kryachokoii 



/Ewald 



(x,y,0) =EE^^(f J )C™oe lmw , (D12) 

j — 1 i.ra 



where Cmo is already introduced by Ea. l45|l and Qj = 
yj{x- s. ]x ) 2 + (y - s jy ) 2 . Applying again the appropri- 
ate (Eg. 145(1 ) representation of the Yl functions, the orig- 
inal equation yields a second order linear inhomogeneous 
ordinary differential equation for fj^, 



f'l + -f' L 



-h 



Af e \2t+l (£ + |m|)! 1 

T 



4tt (£ - |m|)! 2l m l|m| 



£ a^,m)(-l) fc e -^y fc+M 



k=0 



P 



,2A+1 ' 



(D13) 



where 2A = i - 2k - \m\ . Only the mod(€ - \m\ , 2) = 
case has nonzero contribution. The general solutions of 
the homogeneous equation are 



m = 



fe,\ m \(Q)= A m (p?) Hm| +B m (w) |m| , |m|>0 



(D14) 



The solutions of the inhomogeneous equation for both 
m = and \m\ > can be obtained by the constant 
variation method. For example, for m = 



A (q) = 

BoG?) = 
which results 

jEwlad 



i 



1 ( Ya{- V 2 q 2 ) 



log (pq) 



(D15) 



(x,y,K) = 



3=1 

A 



l0g(p 2 Q 2 ) 



(D16) 



by the charge neutrality. Because of the identity (see 
GradshteynM 8.212/1.) 

T(0, x)-^L log (a) = - log (x) (l + A " 



7- 



4?r 
t 



dt 



t 



x > 0, 



(D17) 



the choice of Ao = — v47r together with the charge neu- 
trality ensures the required vanishing potential at the 
infinity (open boundary). 
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As it is shown above in the simple example the ob- 
tained results are the same what are already presented by 
Langridge and co-workersifi or by Ported but of course, 
the followed path is methodologically different. 
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